7-5 Days alive and free of invasive or non-invasive ventilation

This outcome was defined as the number of days alive and free of invasive or non-invasive ventilation from randomisation to 28 days. All patients dying within 28 days will be assigned zero free days.

Authors

Affiliations

James Totterdell

University of Sydney

Rob Mahar

University of Melbourne

Preamble

Load packages

library(tidyverse)
library(labelled)
library(kableExtra)
library(cmdstanr)
library(posterior)
library(bayestestR)
library(bayesplot)
library(matrixStats)
library(plotly)
library(lubridate)
library(ggdist)
library(patchwork)

theme_set(theme_classic(base_size = 10, base_family = "Palatino") +
  theme(panel.grid = element_blank(),
        strip.background = element_blank()))
bayesplot_theme_set(theme_set(theme_classic(base_size = 10, base_family = "Palatino") +
  theme(panel.grid = element_blank(),
        strip.background = element_blank())))

color_scheme_set("red")
options(digits = 4)

Prepare dataset

devtools::load_all()
all_dat <- read_all_no_daily()
all_dat_daily <- read_all_daily()

acs_itt_dat <- all_dat %>% 
  filter_acs_itt() %>%
  transmute_model_cols_grp_aus_nz()
acs_itt_nona_dat <- acs_itt_dat %>% 
  filter(!is.na(out_dafv))

# Concurrent enrolments for C2
acs_itt_concurc2_dat <- acs_itt_dat %>%
  filter_concurrent_intermediate()
acs_itt_concurc2_nona_dat <- acs_itt_concurc2_dat %>% 
  filter(!is.na(out_dafv))

# Concurrent enrolments for C3
acs_itt_concurc3_dat <- acs_itt_dat %>%
  filter_concurrent_std_aspirin()
acs_itt_concurc3_nona_dat <- acs_itt_concurc3_dat %>% 
  filter(!is.na(out_dafv))

# Concurrent enrolments for C4
acs_itt_concurc4_dat <- acs_itt_dat %>%
  filter_concurrent_therapeutic()
acs_itt_concurc4_nona_dat <- acs_itt_concurc4_dat %>% 
  filter(!is.na(out_dafv))

Load models

ordmod0 <- cmdstan_model(
  "stan/ordinal/logistic_cumulative.stan") # No epoch or site
ordmod_epoch <- cmdstan_model(
  "stan/ordinal/logistic_cumulative_epoch.stan") # Epoch only
ordmod <- cmdstan_model(
  "stan/ordinal/logistic_cumulative_epoch_site.stan") # Full model

Code

make_summary_table <- function(dat, format = "html") {
  tdat <- dat %>%
    group_by(CAssignment = factor(
      CAssignment, 
      levels = c("C1", "C2", "C3", "C4"),
      labels = intervention_labels2()$CAssignment[-1])) %>%
    summarise(
      Patients = n(),
      Known = sum(!is.na(out_dafv)),
      Deaths = sprintf("%i (%.0f%%)", sum(D28_death, na.rm = TRUE), 100 * mean(D28_death, na.rm = TRUE)),
      `Any ventilation` = sprintf("%i (%.0f%%)", 
                                  sum(D28_OutcomeDaysFreeOfVentilation < 28, na.rm = TRUE),
                                  100 * mean(D28_OutcomeDaysFreeOfVentilation < 28, na.rm = TRUE)),
      `DAFV, Median (Q1, Q3)` = sprintf(
        "%.0f (%.0f, %.0f)", 
        median(out_dafv, na.rm = T), 
        quantile(out_dafv, 0.25, na.rm = TRUE), 
        quantile(out_dafv, 0.75, na.rm = TRUE))
    ) %>%
    bind_rows(
      dat %>%
    group_by(CAssignment = "Overall") %>%
    summarise(
      Patients = n(),
      Known = sum(!is.na(out_dafv)),
      Deaths = sprintf("%i (%.0f%%)", sum(D28_death, na.rm = TRUE), 100 * mean(D28_death, na.rm = TRUE)),
      `Any ventilation` = sprintf("%i (%.0f%%)", 
                              sum(D28_OutcomeDaysFreeOfVentilation < 28, na.rm = TRUE),
                              100 * mean(D28_OutcomeDaysFreeOfVentilation < 28, na.rm = TRUE)),
      `DAFV, Median (Q1, Q3)` = sprintf(
        "%.0f (%.0f, %.0f)", 
        median(out_dafv, na.rm = T), 
        quantile(out_dafv, 0.25, na.rm = TRUE), 
        quantile(out_dafv, 0.75, na.rm = TRUE))
    )
    ) %>%
    rename(`Anticoagulation\nintervention` = CAssignment)
  kable(
    tdat,
    format = format,
    align = "lrrrrr",
    booktabs = TRUE,
    linesep = ""
  ) %>%
    kable_styling(
      font_size = 9,
      latex_options = "HOLD_position"
    ) %>%
    row_spec(nrow(tdat), bold = T)
}


make_primary_model_data <- function(
    dat, 
    vars = NULL,
    beta_sd_var = NULL,
    ctr = contr.orthonorm,
    p_mult = 2) {
  
  X <- make_X_design(dat, vars, ctr)
  attX <- attr(X, "contrasts")
  X <- X[, -1]
  attr(X, "contrasts") <- attX
  nXtrt <- sum(grepl("rand", colnames(X)))
  
  beta_sd <- c(rep(1, nXtrt), beta_sd_var)
  
  epoch <- dat$epoch
  M_epoch <- max(dat$epoch)
  region <- dat$ctry_num
  M_region <- max(region)
  site <- dat$site_num
  M_site <- max(site)
  region_by_site <- dat %>% 
    dplyr::count(ctry_num, site_num) %>% 
    pull(ctry_num)

  y_raw <- ordered(dat$out_dafv)
  y_mod <- as.integer(y_raw)
  N <- dim(X)[1]
  K <- dim(X)[2]
  unique_y <- length(unique(y_mod))
  list(N = N, K = K, J = unique_y, X = X, y = y_mod, y_raw = y_raw,
       M_region = M_region,
       M_site = M_site, site = site,
       M_epoch = M_epoch, epoch = epoch,
       region_by_site = region_by_site,
       p_par = p_mult * rep(1 / unique_y, unique_y),
       beta_sd = beta_sd)  
}

Outcome Derivation

This outcome was defined as the number of days alive and free of positive pressure ventilation (invasive or non-invasive) to 28 days. All patients dying within 28 days will be assigned zero free days.

The outcome is calculated for a patient as:

• missing, if day 28 mortality is unknown (D28_PatientStatus) or if the patient was known to have been alive at day 28 but number of days requiring ventilation (D28_OutcomeDaysFreeOfVentilation) is unknown
• 0 if they died by day 28 (D28_PatientStatus)
• min(28, D28_OutcomeDaysFreeOfVentilation) otherwise

Descriptive Analyses

Overall

The overall distribution of days alive and free of ventilation (DAFV) are reported in Figure 1 for all participants in the ACS-ITT set.

There were few participants who required any ventilation during hospital who did not die (DAFV of 0),

Code

pdat <- acs_itt_dat %>%
  dplyr::count(dafv = addNA(ordered(out_dafv, levels = 0:28)), .drop = F) %>%
  mutate(p = n / sum(n))
ggplot(pdat, aes(dafv, p)) +
  geom_col() +
  labs(
    x = "Days alive and free of ventilation", 
    y = "Proportion"
  )

(a) ACS-ITT

Figure 1: Distribution of days alive and free of hospital.

Code

save_tex_table(
  make_summary_table(acs_itt_dat, "latex"), 
  file.path("outcomes", "secondary", "7-5-anticoagulation-summary"))
make_summary_table(acs_itt_dat)

Anticoagulation intervention	Patients	Known	Deaths	Any ventilation	DAFV, Median (Q1, Q3)
Low-dose	610	596	19 (3%)	34 (6%)	28 (28, 28)
Intermediate-dose	613	603	15 (2%)	23 (4%)	28 (28, 28)
Low-dose with aspirin	283	281	10 (4%)	18 (6%)	28 (28, 28)
Therapeutic-dose	50	50	6 (12%)	7 (14%)	28 (28, 28)
Overall	1556	1530	50 (3%)	82 (5%)	28 (28, 28)

Summaty of DAFV.

Anti-coagulation

Code

pdat <- acs_itt_nona_dat %>%
  dplyr::count(
    CAssignment = factor(CAssignment, labels = str_replace(intervention_labels()$CAssignment[-1], "<br>", "\n")),
    dafv = ordered(out_dafv, levels = 0:28), 
    .drop = F) %>%
  group_by(CAssignment) %>%
  mutate(p = n / sum(n)) %>%
  mutate(cp = cumsum(p)) %>%
  ungroup()
p <- ggplot(pdat, aes(CAssignment, p)) +
  geom_col(aes(fill = fct_rev(dafv))) +
  labs(x = "Anticoagulation intervention", y = "Cumulative proportion") +
  scale_fill_viridis_d("Days alive and\nfree of ventilation", option = "A", direction = -1) +
  guides(fill = guide_legend(reverse = TRUE)) +
  theme(legend.key.size = unit(0.75, "lines"))
pth <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(pth, "outcome-7-5-descriptive.pdf"), p, height = 3, width = 6)
p

Figure 2: Distribution of days alive and free of hospital by anti-coagulation intervention.

Age

DAFV by age

pdat <- acs_itt_nona_dat %>%
  dplyr::count(
    agegrp = cut(AgeAtEntry, c(18, seq(25, 75, 5), 100), include.lowest = T, right = F),
    dafv = fct_rev(ordered(out_dafv, levels = 0:28)), 
    .drop = F) %>%
  group_by(agegrp) %>%
  mutate(p = n / sum(n))
pdat2 <- pdat %>%
  group_by(agegrp) %>%
  summarise(n = sum(n))
p1 <- ggplot(pdat2, aes(agegrp, n)) +
  geom_col(colour = "grey40", fill = "grey40") +
  geom_vline(xintercept = 60, linetype = 2) +
  labs(y = "Number of\nparticipants",
       x = "Age at entry") +
  geom_vline(xintercept = 8.5, linetype = 2) +
  theme(axis.text.x = element_text(hjust = 1, angle = 45))
p2 <- ggplot(pdat, aes(agegrp, p)) +
  geom_col(aes(fill = dafv)) +
  labs(x = "Age", y = "Cumulative\nproportion") +
  scale_fill_viridis_d("DAFV", option = "A", direction = -1) +
  guides(fill = guide_legend(reverse = TRUE)) +
  geom_vline(xintercept = 8.5, linetype = 2) +
  theme(axis.text.x = element_text(hjust = 1, angle = 45)) +
  theme(legend.key.size = unit(0.25, "lines"))
p <- p1 | p2
path <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(path, "7-5-age.pdf"), p, height = 2.5, width = 6)
p

Figure 3: Distribution of days alive and free of ventilation by age group

Country

DAFV by country

pdat <- acs_itt_nona_dat %>%
  dplyr::count(
    country,
    dafv = ordered(out_dafv, levels = 0:28), 
    .drop = F) %>%
  group_by(country) %>%
  mutate(p = n / sum(n)) %>%
  mutate(cp = cumsum(p)) %>%
  ungroup()
pdat2 <- pdat %>%
  group_by(country) %>%
  summarise(n = sum(n))

p1 <- ggplot(pdat2, aes(country, n)) +
  geom_col() +
    labs(
      y = "Number of\nparticipants", 
      x = "Country of enrolment")

p2 <- ggplot(pdat, aes(country, p)) +
  geom_col(aes(fill = fct_rev(dafv))) +
  labs(x = "Anticoagulation intervention", y = "Cumulative\nproportion") +
  scale_fill_viridis_d("DAFV", option = "A", direction = -1) +
  guides(fill = guide_legend(reverse = TRUE)) +
  theme(legend.key.size = unit(0.25, "lines"))
p <- p1 | p2
path <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(path, "7-5-country.pdf"), p, height = 2.5, width = 6)
p

Figure 4: Distribution of days alive and free of ventilation by country.

Site

DAFV by site

pdat <- all_dat %>%
  filter_acs_itt() %>%
  filter(!is.na(out_dafv)) %>%
  dplyr::count(
    Country = factor(PT_CountryName, levels = c("India", "Australia", "Nepal", "New Zealand"),
                     labels = c("India", "Australia", "Nepal", "New\nZealand")),
    Site = fct_infreq(Location),
    dafv = ordered(out_dafv, levels = 0:28)) %>%
  complete(dafv = ordered(0:28), nesting(Country, Site), fill = list(n = 0)) %>%
  group_by(Country, Site) %>%
  mutate(p = n / sum(n)) %>%
  mutate(cp = cumsum(p)) %>%
  ungroup() %>%
  mutate(
    Country = droplevels(Country),
    Site = droplevels(Site)
  )
pdat2 <- pdat %>%
  group_by(Country, Site) %>%
  summarise(n = sum(n)) %>%
  ungroup()
p1 <- ggplot(pdat2, aes(Site, n)) +
  facet_grid( ~ Country, scales = "free_x", space = "free_x") +
  geom_col() +
    labs(
      y = "Number of\nparticipants", 
      x = "") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        panel.border = element_rect(fill = NA))
p2 <- ggplot(pdat, aes(Site, p)) +
  facet_grid( ~ Country, scales = "free_x", space = "free_x") +
  geom_col(aes(fill = fct_rev(dafv))) +
  labs(x = "Anticoagulation intervention", y = "Cumulative\nproportion") +
  scale_fill_viridis_d("DAFV", option = "A", direction = -1) +
  guides(fill = guide_legend(reverse = TRUE, ncol = 1)) +
  theme(axis.text.x = element_text(hjust = 1, angle = 45)) +
  theme(legend.key.size = unit(0.25, "lines"))
p <- p1 / p2 +
  plot_layout(guides = 'collect')
path <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(path, "7-5-country-site.pdf"), p, height = 4, width = 6.25)
p

Figure 5: Distribution of days alive and free of ventilation by study site

Calendar Time

DAFV by calendar date

pdat <- acs_itt_nona_dat %>%
  dplyr::count(
    yr = year(RandDate), mth = month(RandDate),
    dafv = ordered(out_dafv, levels = 0:28), 
    .drop = F) %>%
  group_by(yr, mth) %>%
  mutate(p = n / sum(n)) %>%
  mutate(cp = cumsum(p)) %>%
  ungroup()
p1 <- pdat %>%
  group_by(yr, mth) %>%
  summarise(n = sum(n)) %>%
  ggplot(., aes(mth, n))  +
  facet_grid( ~ yr, drop = T, scales = "free_x", space = "free") +
    geom_col() +
    labs(
      y = "Number of\nparticipants", 
      x = "Calendar date (month of year)") +
  scale_x_continuous(breaks = 1:12)
p2 <- ggplot(pdat, aes(mth, p)) +
  facet_grid( ~ yr, drop = T, scales = "free_x", space = "free") +
  geom_col(aes(fill = fct_rev(dafv))) +
  labs(x = "Calendar date (month of year)", y = "Cumulative\nproportion") +
  scale_fill_viridis_d("DAFV", option = "A", direction = -1) +
  guides(fill = guide_legend(reverse = TRUE)) +
  theme(legend.key.size = unit(0.25, "lines")) +
  scale_x_continuous(breaks = 1:12)
p <- p1 | p2
path <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(path, "7-5-calendar-time.pdf"), p, height = 2, width = 6)
p

Figure 6: Distribution of days alive and free of ventilation by calendar time.

Sample Cumulative Logits

Plot sample cumulative logits

trt_counts <- acs_itt_nona_dat %>%
  dplyr::count(CAssignment, out_dafv) %>%
  complete(CAssignment, out_dafv, fill = list(n = 0)) %>%
  group_by(CAssignment) %>%
  mutate(p = n / sum(n))
trt_logit <- trt_counts %>% 
  group_by(CAssignment) %>% 
  mutate(clogit = logit(cumsum(p))) %>%
  group_by(out_dafv) %>%
  mutate(rel_clogit = clogit - mean(clogit)) %>%
  filter(out_dafv != 28)

ggplot(trt_logit, aes(CAssignment, rel_clogit)) +
  facet_wrap( ~ out_dafv) +
  geom_point() +
  labs(y = "Relative (to mean) sample cumulative logit")

Inspect sample cumulative logits.

Plot sample cumulative logits

ggplot(trt_logit, aes(out_dafv, rel_clogit)) +
  facet_wrap( ~ CAssignment) +
  geom_point() +
  labs(y = "Relative (to mean) sample cumulative logit")

Inspect sample cumulative logits.

Modelling

The primary analysis uses the ACS-ITT set, so all participants have been randomised to the anticoagulation domain and some participants may have been randomised to the antiviral domain.

Primary Model

Fit primary model

fit_primary_model <- function(dat, ctr = contr.orthonorm, save = FALSE) {
  mdat <- make_primary_model_data(
    dat, c("inelgc3", "agegte60", "ctry"), 
    c(10, 2.5, 1, 1), 
    ctr)
  snk <- capture.output(
  mfit <- ordmod$sample(
    data = mdat,
    seed = 3510392,
    chains = 8,
    parallel_chains = 8,
    iter_warmup = 1000,
    iter_sampling = 2500,
    refresh = 0,
    adapt_delta = 0.98
  )) 
  if (save) {
    mfit$save_object(file.path("outputs", "models", "secondary", "7-5.rds"))    
  }
  mdrws <- as_draws_rvars(mfit$draws())
  # Add names
  epoch_map <- dat %>% dplyr::count(epoch, epoch_lab)
  site_map <- dat %>% dplyr::count(site_num, site)
  names(mdrws$beta) <- colnames(mdat$X)
  names(mdrws$gamma_epoch) <- epoch_map$epoch_lab
  names(mdrws$gamma_site) <- site_map$site
  # Convert to treatment log odds ratios
  mdrws$Acon <- attr(mdat$X, "contrasts")$randA %**% mdrws$beta[grepl("randA[0-9]", names(mdrws$beta))]
  mdrws$Ccon <- attr(mdat$X, "contrasts")$randC %**% mdrws$beta[grepl("randC[0-9]", names(mdrws$beta))]
  mdrws$trtA <- mdrws$Acon[-1] - mdrws$Acon[1]
  mdrws$trtC <- mdrws$Ccon[-1] - mdrws$Ccon[1]
  mdrws$compare <- c(
    "Intermediate vs low" = exp(mdrws$trtC[1]),
    "Low with aspirin vs low" = exp(mdrws$trtC[2]),
    "Therapeutic vs low" = exp(mdrws$trtC[3]),
    "Intermediate vs low with aspirin" = exp(mdrws$trtC[1] - mdrws$trtC[2]),
    "Intermediate vs therapeutic" = exp(mdrws$trtC[1] - mdrws$trtC[3]),
    "Low with aspirin vs therapeutic" = exp(mdrws$trtC[2] - mdrws$trtC[3])
  )
  mdrws$OR <- c(
    setNames(exp(mdrws$trtC), c("Intermediate", "Low with aspirin", "Therapeutic")),
    "Ineligible aspirin" = exp(mdrws$beta[which(grepl("inelg", colnames(mdat$X)))]),
    "Age \u2265 60" = exp(mdrws$beta[which(grepl("age", colnames(mdat$X)))]),
    setNames(exp(mdrws$beta[which(grepl("ctry", colnames(mdat$X)))]), c("AU/NZ", "NP"))
  )
  return(list(dat = mdat, fit = mfit, drws = mdrws))
}
res <- fit_primary_model(acs_itt_nona_dat, save = TRUE)

Odds ratio summary table

save_tex_table(
  odds_ratio_summary_table_rev(res$drws$OR, "latex"),
  "outcomes/secondary/7-5-primary-model-acs-itt-summary-table")
odds_ratio_summary_table_rev(res$drws$OR)

Parameter	Median	95% CrI	Mean (SD)	Pr(OR > 1)
Intermediate	1.40	(0.81, 2.45)	1.46 (0.42)	0.88
Low with aspirin	1.16	(0.62, 2.19)	1.22 (0.41)	0.68
Therapeutic	0.41	(0.15, 1.15)	0.47 (0.27)	0.05
Ineligible aspirin	0.29	(0.09, 1.08)	0.37 (0.28)	0.03
Age ≥ 60	0.50	(0.31, 0.83)	0.52 (0.13)	0.00
AU/NZ	1.10	(0.27, 4.92)	1.47 (1.28)	0.55
NP	0.56	(0.15, 2.32)	0.74 (0.68)	0.20

Posterior summaries for model parameters (fixed-effects).

Odds ratio summary for epoch and site

p <- plot_epoch_site_terms(
  exp(res$drws$gamma_epoch),
  exp(res$drws$gamma_site),
  factor(res$dat$region_by_site, 
         labels = c("India", "Australia\nNew Zealand", "Nepal"))
)
pth <- file.path("outputs", "figures", "outcomes", "secondary", "7-5-primary-model-epoch-site-terms.pdf")
ggsave(pth, p, width = 6, height = 4.5)
p

Summary of epoch and site posterior odds ratios.

Code

p <- plot_or_densities(res$drws$compare)
p

Posterior densities for treatment comparisons.

Diagnostics

Code

res$fit$summary(c("alpha", "beta", "gamma_site", "gamma_epoch", "tau_site", "tau_epoch")) %>%
  print(n = Inf)

# A tibble: 68 × 10
   variable            mean   median    sd   mad      q5     q95  rhat ess_bulk
   <chr>              <dbl>    <dbl> <dbl> <dbl>   <dbl>   <dbl> <dbl>    <dbl>
 1 alpha[1]        -4.37    -4.36    0.704 0.682 -5.55   -3.23    1.00    2883.
 2 alpha[2]        -4.32    -4.32    0.704 0.684 -5.50   -3.18    1.00    2869.
 3 alpha[3]        -4.30    -4.29    0.704 0.682 -5.48   -3.16    1.00    2858.
 4 alpha[4]        -4.28    -4.27    0.704 0.682 -5.46   -3.14    1.00    2850.
 5 alpha[5]        -4.26    -4.25    0.703 0.679 -5.43   -3.11    1.00    2842.
 6 alpha[6]        -4.23    -4.23    0.703 0.679 -5.41   -3.10    1.00    2840.
 7 alpha[7]        -4.19    -4.19    0.702 0.678 -5.37   -3.06    1.00    2832.
 8 alpha[8]        -4.17    -4.17    0.702 0.677 -5.35   -3.03    1.00    2823.
 9 alpha[9]        -4.13    -4.13    0.702 0.679 -5.30   -3.00    1.00    2812.
10 alpha[10]       -4.10    -4.09    0.701 0.677 -5.27   -2.97    1.00    2809.
11 alpha[11]       -4.08    -4.07    0.701 0.677 -5.24   -2.95    1.00    2805.
12 alpha[12]       -4.06    -4.05    0.701 0.678 -5.22   -2.93    1.00    2800.
13 alpha[13]       -4.02    -4.01    0.701 0.678 -5.19   -2.89    1.00    2793.
14 alpha[14]       -3.99    -3.98    0.700 0.677 -5.16   -2.86    1.00    2782.
15 alpha[15]       -3.94    -3.93    0.700 0.676 -5.11   -2.81    1.00    2778.
16 alpha[16]       -3.86    -3.85    0.699 0.674 -5.03   -2.74    1.00    2780.
17 alpha[17]       -3.81    -3.80    0.699 0.673 -4.97   -2.68    1.00    2771.
18 beta[1]         -0.753   -0.758   0.420 0.426 -1.44   -0.0600  1.00   18460.
19 beta[2]          0.118    0.117   0.226 0.228 -0.251   0.493   1.00   20145.
20 beta[3]          0.555    0.553   0.279 0.281  0.0984  1.02    1.00   18259.
21 beta[4]         -0.121   -0.115   0.694 0.688 -1.28    1.00    1.00   16789.
22 beta[5]          0.893    0.882   0.530 0.526  0.0392  1.79    1.00   23035.
23 beta[6]         -1.21    -1.23    0.633 0.630 -2.22   -0.134   1.00   24360.
24 beta[7]         -0.691   -0.692   0.255 0.254 -1.10   -0.269   1.00   19969.
25 beta[8]          0.112    0.0971  0.732 0.729 -1.07    1.35    1.00   10389.
26 beta[9]         -0.560   -0.578   0.692 0.676 -1.67    0.594   1.00    6683.
27 gamma_site[1]    2.07     1.82    1.50  1.27   0.139   4.83    1.00    8903.
28 gamma_site[2]    2.25     2.00    1.51  1.27   0.331   4.96    1.00    7889.
29 gamma_site[3]    1.48     1.25    1.54  1.31  -0.534   4.26    1.00   10176.
30 gamma_site[4]   -1.06    -1.08    0.601 0.572 -2.02   -0.0458  1.00    4793.
31 gamma_site[5]   -0.0208  -0.0425  0.596 0.565 -0.956   0.996   1.00    4206.
32 gamma_site[6]   -1.61    -1.62    0.626 0.610 -2.63   -0.582   1.00    5828.
33 gamma_site[7]    2.27     2.02    1.52  1.27   0.346   5.06    1.00    8546.
34 gamma_site[8]   -0.440   -0.450   0.692 0.666 -1.56    0.707   1.00    6550.
35 gamma_site[9]    1.47     1.40    0.848 0.808  0.221   2.97    1.00    7327.
36 gamma_site[10]   1.98     1.74    1.50  1.28   0.0214  4.74    1.00    8735.
37 gamma_site[11]  -0.387   -0.400   0.657 0.628 -1.43    0.716   1.00    5897.
38 gamma_site[12]  -0.569   -0.435   1.03  0.890 -2.41    0.933   1.00   16270.
39 gamma_site[13]   0.627    0.406   1.18  0.866 -0.898   2.78    1.00   19907.
40 gamma_site[14]  -1.05    -0.906   1.04  1.05  -2.94    0.361   1.00   11400.
41 gamma_site[15]   0.534    0.333   1.21  0.884 -1.09    2.70    1.00   19995.
42 gamma_site[16]  -0.989   -0.875   0.988 1.02  -2.78    0.369   1.00   11863.
43 gamma_site[17]  -0.976   -0.894   0.887 0.905 -2.57    0.253   1.00   11363.
44 gamma_site[18]   0.425    0.232   1.21  0.900 -1.23    2.61    1.00   21990.
45 gamma_site[19]   0.336    0.166   1.22  0.906 -1.39    2.49    1.00   21820.
46 gamma_site[20]  -0.587   -0.459   1.03  0.892 -2.44    0.902   1.00   16073.
47 gamma_site[21]   0.438    0.246   1.22  0.901 -1.24    2.63    1.00   22541.
48 gamma_site[22]   0.502    0.294   1.23  0.903 -1.15    2.71    1.00   21711.
49 gamma_site[23]   1.15     0.938   1.15  1.00  -0.233   3.29    1.00   14574.
50 gamma_site[24]   0.625    0.409   1.18  0.872 -0.924   2.80    1.00   18145.
51 gamma_site[25]  -0.270   -0.126   0.579 0.359 -1.39    0.454   1.00   11367.
52 gamma_site[26]  -0.00143  0.00867 0.557 0.335 -0.965   0.859   1.00   15002.
53 gamma_epoch[1]   0        0       0     0      0       0      NA         NA 
54 gamma_epoch[2]  -0.229   -0.169   0.402 0.332 -0.969   0.338   1.00   10447.
55 gamma_epoch[3]  -0.333   -0.268   0.478 0.427 -1.20    0.345   1.00    8933.
56 gamma_epoch[4]  -0.266   -0.224   0.464 0.415 -1.08    0.445   1.00    9924.
57 gamma_epoch[5]  -0.545   -0.474   0.561 0.551 -1.58    0.231   1.00    7721.
58 gamma_epoch[6]  -0.373   -0.325   0.554 0.507 -1.33    0.469   1.00   10405.
59 gamma_epoch[7]  -0.480   -0.430   0.569 0.550 -1.48    0.357   1.00    8708.
60 gamma_epoch[8]  -0.567   -0.515   0.577 0.564 -1.58    0.283   1.00    8024.
61 gamma_epoch[9]  -0.795   -0.751   0.608 0.617 -1.86    0.0821  1.00    6530.
62 gamma_epoch[10] -1.08    -1.05    0.680 0.698 -2.26   -0.0280  1.00    5706.
63 gamma_epoch[11] -0.595   -0.550   0.616 0.616 -1.66    0.328   1.00    7632.
64 gamma_epoch[12] -0.307   -0.281   0.645 0.586 -1.40    0.741   1.00   10145.
65 tau_site[1]      1.83     1.68    0.733 0.603  0.937   3.22    1.00    5368.
66 tau_site[2]      1.20     1.12    0.680 0.616  0.231   2.42    1.00    6674.
67 tau_site[3]      0.647    0.483   0.619 0.454  0.0438  1.80    1.00   11241.
68 tau_epoch        0.439    0.412   0.237 0.224  0.0916  0.863   1.00    5973.
# … with 1 more variable: ess_tail <dbl>

Code

res$fit$diagnostic_summary()

$num_divergent
[1] 0 0 0 0 0 0 0 0

$num_max_treedepth
[1] 0 0 0 0 0 0 0 0

$ebfmi
[1] 0.7775 0.7944 0.7418 0.8248 0.7489 0.7897 0.8120 0.7979

Trace plots

Code

mcmc_trace(res$drws["beta"])

Code

mcmc_trace(res$drws["alpha"])

Code

mcmc_trace(res$drws["gamma_site"])

Code

mcmc_trace(res$drws["gamma_epoch"])

Posterior Predictive

Code

y_raw <- as.integer(levels(res$dat$y_raw))
y_ppc <- res$drws$y_ppc
y_ppc_raw <- rfun(\(x) y_raw[x])(y_ppc)
ppc_dat <- bind_cols(acs_itt_nona_dat, tibble(y_ppc = y_ppc_raw)) %>%
  mutate(CAssignment = factor(
    CAssignment, 
    levels = names(intervention_labels_short_break()$CAssignment),
    labels = intervention_labels_short_break()$CAssignment))

grp_ppc <- function(grp, d = 0:28) {
  ppc_dat %>%
    group_by(grp = {{grp}}) %>%
    summarise(
      y_lteq = map_dbl(d, ~ mean(out_dafv <= .x)),
      ypp_lteq = map(d, ~ rvar_mean(y_ppc <= .x)),
      y_eq = map_dbl(d, ~ mean(out_dafv == .x)),
      ypp_eq = map(d, ~ rvar_mean(y_ppc == .x))
    ) %>%
    unnest(c(ypp_lteq, ypp_eq)) %>%
    mutate(days = d) %>%
    ungroup() %>%
    pivot_longer(
      y_lteq:ypp_eq, 
      names_to = c("response", "event"),
      names_sep = "_",
      values_to = "posterior")
}

plot_grp_ppc <- function(dat, lab = "", xlab = "Probability") {
  ggplot(dat %>% 
           filter(response == "ypp", event == "eq"), 
         aes(y = days)) +
    facet_wrap( ~ grp, nrow = 1, scales = "free_x") +
    stat_slabinterval(aes(xdist = posterior), fatten_point = 1)  +
    geom_point(data = dat %>% filter(response == "y", event == "eq"), 
               aes(y = days, x = mean(posterior)),
               colour = "red",
               shape = 23) +
    scale_x_continuous(
      xlab, breaks = c(0, 0.3, 0.6),
      sec.axis = sec_axis(
        ~ . , name = lab, breaks = NULL, labels = NULL)) +
    scale_y_continuous("Days") +
    theme(strip.text = element_text(size = rel(0.7)),
          axis.title.x = element_text(size = rel(0.7)),
          axis.text.x = element_text(size = rel(0.65)),
          axis.title.y = element_text(size = rel(0.75)),
          axis.title.x.bottom = element_blank())
}

pp_C <- grp_ppc(CAssignment)
pp_ctry <- grp_ppc(country)
pp_epoch <- grp_ppc(epoch)
pp_site <- grp_ppc(site) %>%
  left_join(ppc_dat %>% dplyr::count(site, country), 
            by = c("grp" = "site"))

p1 <- plot_grp_ppc(pp_C, "Anticoagulation", "")
p2 <- plot_grp_ppc(pp_ctry, "Country", "") 
p3 <- plot_grp_ppc(pp_epoch, "Epoch", "")
p4 <- plot_grp_ppc(
  pp_site %>% filter(country == "IN"), "Sites India", "")
p5 <- plot_grp_ppc(
  pp_site %>% filter(country == "AU"), "Sites Australia", "")
p6 <- plot_grp_ppc(
  pp_site %>% filter(country == "NP"), "Sites Nepal", "")
p7 <- plot_grp_ppc(
  pp_site %>% filter(country == "NP"), "Sites New Zealand", "")

g1 <- (p1 | p2) / p3
g2 <- p4 / p5 / (p6 | p7)

pth1 <- file.path(
  "outputs", "figures", "outcomes", "secondary",
  "7-5-primary-model-acs-itt-ppc1.pdf")
pth2 <- file.path(
  "outputs", "figures", "outcomes", "secondary",
  "7-5-primary-model-acs-itt-ppc2.pdf")

ggsave(pth1, g1, width = 6, height = 5, device = cairo_pdf)
ggsave(pth2, g2, width = 6, height = 6, device = cairo_pdf)

g1

Code

g2

Sensitivity: Concurrent controls

Three separate analyses are conducted based on concurrent randomisations as was done for the primary outcome. For these models, the epoch term has been removed.

Intermediate dose

Descriptive summary table

save_tex_table(
  make_summary_table(acs_itt_concurc2_dat, "latex"), 
  file.path("outcomes", "secondary", "7-5-anticoagulation-concurrent-intermediate-summary"))
make_summary_table(acs_itt_concurc2_dat)

Anticoagulation intervention	Patients	Known	Deaths	Any ventilation	DAFV, Median (Q1, Q3)
Low-dose	610	596	19 (3%)	34 (6%)	28 (28, 28)
Intermediate-dose	613	603	15 (2%)	23 (4%)	28 (28, 28)
Overall	1223	1199	34 (3%)	57 (5%)	28 (28, 28)

Fit model

ctr <- contr.orthonorm
X <- model.matrix( 
  ~ CAssignment + agegte60 + ctry, 
  data = acs_itt_concurc2_nona_dat, 
  contrasts.arg = list(CAssignment = ctr))[, -1]
y_raw <- ordered(acs_itt_concurc2_nona_dat$out_dafv)
y_mod <- as.integer(y_raw)
N <- dim(X)[1]
K <- dim(X)[2]
unique_y <- length(unique(y_mod))
mdat <- list(
  N = N, K = K, J = unique_y, X = X, y = y_mod, y_raw = y_raw,
  p_par = 2 * rep(1 / unique_y, unique_y),
  beta_sd = c(1, 2.5, 1, 1)
)
snk <- capture.output(
  mfit <- ordmod0$sample(
    data = mdat,
    seed = 75136,
    chains = 8,
    parallel_chains = 8,
    iter_warmup = 1000,
    iter_sampling = 2500,
    refresh = 0
))
mdrws <- as_draws_rvars(mfit$draws())
names(mdrws$beta) <- colnames(X)
mdrws$Ccon <- as.vector(ctr(2) %**% mdrws$beta[1])
mdrws$Ctrt <- mdrws$Ccon[-1] - mdrws$Ccon[1]
mdrws$compare <- c("Intermediate vs low" = exp(mdrws$Ctrt[1]))
mdrws$OR <- c("Intermediate" = exp(mdrws$Ctrt), 
              setNames(exp(mdrws$beta[2:4]), c("Age \u2265 60", "Australia/New Zealand", "Nepal")))

Sample distribution

pdat <- acs_itt_concurc2_nona_dat %>%
  dplyr::count(
    CAssignment = factor(
      CAssignment, 
      levels = c("C1", "C2"),
      labels = c("Low", "Intermediate")),
    dafv = ordered(out_dafv, levels = 0:28), 
    .drop = F) %>%
  group_by(CAssignment) %>%
  mutate(p = n / sum(n)) %>%
  mutate(cp = cumsum(p)) %>%
  ungroup()
p <- ggplot(pdat, aes(CAssignment, p)) +
  geom_col(aes(fill = fct_rev(dafv)), colour = NA) +
  labs(x = "Anticoagulation intervention", y = "Cumulative proportion") +
  scale_fill_viridis_d("Days alive and\nfree of ventilation", option = "A", direction = -1) +
  guides(fill = guide_legend(reverse = TRUE, ncol = 3)) +
  theme(legend.key.size = unit(0.75, "lines"))
pth <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(pth, "outcome-7-5-descriptive-concurrent-intermediate.pdf"), p, height = 2.5, width = 4)
p

Distribution of days alive and free of hospital by anti-coagulation intervention in subset.

Posterior contrast

p <- plot_or_densities(mdrws$compare) +
  geom_text(aes(x = 1, y = 1 - 0.1, 
                label = sprintf("Pr(OR > 1) = %.2f", Pr(RV > 1))), 
            hjust = 0, nudge_x = 0.01, family = "Palatino") +
  scale_x_log10(breaks = c(0.8, 0.9, 1, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6)) +
  labs(y = "Comparison")
p

Odds ratio summary table

save_tex_table(
  odds_ratio_summary_table_rev(mdrws$OR, "latex"),
  "outcomes/secondary/7-5-primary-model-acs-itt-concurrent-intermediate-summary-table")
odds_ratio_summary_table_rev(mdrws$OR)

Parameter	Median	95% CrI	Mean (SD)	Pr(OR > 1)
Intermediate	1.51	(0.91, 2.59)	1.57 (0.43)	0.94
Age ≥ 60	0.43	(0.25, 0.74)	0.44 (0.13)	0.00
Australia/New Zealand	0.73	(0.36, 1.62)	0.80 (0.33)	0.21
Nepal	0.92	(0.40, 2.44)	1.05 (0.54)	0.43

Posterior summaries for model parameters.

Diagnostics

Code

mfit$summary()

# A tibble: 40 × 10
   variable       mean   median      sd     mad       q5      q95  rhat ess_bulk
   <chr>         <dbl>    <dbl>   <dbl>   <dbl>    <dbl>    <dbl> <dbl>    <dbl>
 1 lp__       -3.78e+2 -3.78e+2 3.29e+0 3.22e+0 -3.84e+2 -3.73e+2  1.00    8006.
 2 beta_raw[…  2.94e-1  2.92e-1 1.90e-1 1.88e-1 -1.79e-2  6.09e-1  1.00   23243.
 3 beta_raw[… -3.40e-1 -3.42e-1 1.11e-1 1.11e-1 -5.23e-1 -1.58e-1  1.00   15990.
 4 beta_raw[… -3.00e-1 -3.12e-1 3.84e-1 3.84e-1 -9.12e-1  3.55e-1  1.00   20207.
 5 beta_raw[… -6.35e-2 -8.54e-2 4.57e-1 4.47e-1 -7.87e-1  7.28e-1  1.00   21016.
 6 p[1]        1.95e-2  1.92e-2 4.30e-3 4.25e-3  1.31e-2  2.71e-2  1.00   16213.
 7 p[2]        6.52e-4  4.69e-4 6.23e-4 4.59e-4  4.28e-5  1.89e-3  1.00   16237.
 8 p[3]        6.58e-4  4.65e-4 6.39e-4 4.57e-4  4.56e-5  1.94e-3  1.00   18456.
 9 p[4]        6.49e-4  4.63e-4 6.25e-4 4.48e-4  4.29e-5  1.89e-3  1.00   14798.
10 p[5]        6.56e-4  4.75e-4 6.24e-4 4.57e-4  4.78e-5  1.87e-3  1.00   17414.
# … with 30 more rows, and 1 more variable: ess_tail <dbl>

Code

mfit$diagnostic_summary()

$num_divergent
[1] 0 0 0 0 0 0 0 0

$num_max_treedepth
[1] 0 0 0 0 0 0 0 0

$ebfmi
[1] 0.9209 1.0093 1.0153 0.9283 0.9173 1.0053 0.9836 0.9998

Trace plots

Code

mcmc_trace(mdrws["beta"])

Code

mcmc_trace(mdrws["alpha"])

Low-dose with aspirin

Descriptive summary table

save_tex_table(
  make_summary_table(acs_itt_concurc3_dat, "latex"), 
  file.path("outcomes", "secondary", "7-5-anticoagulation-concurrent-stdaspirin-summary"))
make_summary_table(acs_itt_concurc3_dat)

Anticoagulation intervention	Patients	Known	Deaths	Any ventilation	DAFV, Median (Q1, Q3)
Low-dose	299	291	11 (4%)	20 (7%)	28 (28, 28)
Intermediate-dose	298	293	12 (4%)	15 (5%)	28 (28, 28)
Low-dose with aspirin	283	281	10 (4%)	18 (6%)	28 (28, 28)
Overall	880	865	33 (4%)	53 (6%)	28 (28, 28)

Fit model

ctr <- contr.orthonorm
X <- model.matrix( 
  ~ CAssignment + agegte60 + ctry, 
  data = acs_itt_concurc3_nona_dat, 
  contrasts.arg = list(CAssignment = ctr))[, -1]
y_raw <- ordered(acs_itt_concurc3_nona_dat$out_dafv)
y_mod <- as.integer(y_raw)
N <- dim(X)[1]
K <- dim(X)[2]
unique_y <- length(unique(y_mod))
mdat <- list(
  N = N, K = K, J = unique_y, X = X, y = y_mod, y_raw = y_raw,
  p_par = 2 * rep(1 / unique_y, unique_y),
  beta_sd = c(1, 1, 2.5, 1)
)
snk <- capture.output(
  mfit <- ordmod0$sample(
    data = mdat,
    seed = 135356,
    chains = 8,
    parallel_chains = 8,
    iter_warmup = 1000,
    iter_sampling = 2500,
    refresh = 0
))
mdrws <- as_draws_rvars(mfit$draws())
names(mdrws$beta) <- colnames(X)
mdrws$Ccon <- as.vector(ctr(3) %**% mdrws$beta[1:2])
mdrws$Ctrt <- mdrws$Ccon[-1] - mdrws$Ccon[1]
mdrws$compare <- c("Intermediate vs low" = exp(mdrws$Ctrt[1]),
                   "Low with aspirin vs low" = exp(mdrws$Ctrt[2]),
                   "Intermediate vs low with aspirin" = exp(mdrws$Ctrt[1] - mdrws$Ctrt[2]))
mdrws$OR <- c("Intermediate" = exp(mdrws$Ctrt[1]),
              "Low with aspirin" = exp(mdrws$Ctrt[2]),
              setNames(exp(mdrws$beta[3:4]), c("Age \u2265 60", "Australia/New Zealand")))

Sample distribution

pdat <- acs_itt_concurc3_nona_dat %>%
  dplyr::count(
    CAssignment = factor(
      CAssignment, 
      levels = c("C1", "C2", "C3"),
      labels = c("Low", "Intermediate", "Low\nwith aspirin")),
    dafv = ordered(out_dafv, levels = 0:28), 
    .drop = F) %>%
  group_by(CAssignment) %>%
  mutate(p = n / sum(n)) %>%
  mutate(cp = cumsum(p)) %>%
  ungroup()
p <- ggplot(pdat, aes(CAssignment, p)) +
  geom_col(aes(fill = fct_rev(dafv))) +
  labs(x = "Anticoagulation intervention", y = "Cumulative proportion") +
  scale_fill_viridis_d("Days alive and\nfree of ventilation", option = "A", direction = -1) +
  guides(fill = guide_legend(reverse = TRUE, ncol = 3)) +
  theme(legend.key.size = unit(0.75, "lines"))
p

Distribution of days alive and free of ventilation by anti-coagulation intervention in subset.

Sample distribution

pth <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(pth, "outcome-7-5-descriptive-concurrent-stdaspirin.pdf"), p, height = 2.5, width = 6)

Posterior contrast

plot_or_densities(mdrws$compare)

Odds ratio summary table

save_tex_table(
  odds_ratio_summary_table_rev(mdrws$OR, "latex"),
  "outcomes/secondary/7-5-primary-model-acs-itt-concurrent-stdaspirin-summary-table")
odds_ratio_summary_table_rev(mdrws$OR)

Parameter	Median	95% CrI	Mean (SD)	Pr(OR > 1)
Intermediate	1.32	(0.68, 2.61)	1.40 (0.50)	0.79
Low with aspirin	1.13	(0.60, 2.19)	1.20 (0.41)	0.65
Age ≥ 60	0.51	(0.29, 0.89)	0.53 (0.16)	0.01
Australia/New Zealand	3.15	(0.90, 14.44)	4.28 (3.96)	0.96

Posterior summaries for model parameters.

Diagnostics

Code

mfit$summary()

# A tibble: 38 × 10
   variable       mean   median      sd     mad       q5      q95  rhat ess_bulk
   <chr>         <dbl>    <dbl>   <dbl>   <dbl>    <dbl>    <dbl> <dbl>    <dbl>
 1 lp__       -3.30e+2 -3.30e+2 3.19    3.10e+0 -3.36e+2 -3.25e+2  1.00    7709.
 2 beta_raw[… -1.05e-1 -1.03e-1 0.247   2.43e-1 -5.18e-1  3.00e-1  1.00   21184.
 3 beta_raw[… -1.65e-1 -1.67e-1 0.235   2.36e-1 -5.53e-1  2.25e-1  1.00   23187.
 4 beta_raw[… -2.70e-1 -2.71e-1 0.115   1.15e-1 -4.57e-1 -8.10e-2  1.00   15091.
 5 beta_raw[…  1.18e+0  1.15e+0 0.710   7.06e-1  8.24e-2  2.42e+0  1.00   19514.
 6 p[1]        3.19e-2  3.14e-2 0.00662 6.52e-3  2.18e-2  4.34e-2  1.00   16106.
 7 p[2]        2.02e-3  1.72e-3 0.00139 1.19e-3  3.96e-4  4.72e-3  1.00   20516.
 8 p[3]        1.08e-3  7.76e-4 0.00103 7.48e-4  7.35e-5  3.15e-3  1.00   17492.
 9 p[4]        1.08e-3  7.71e-4 0.00103 7.55e-4  7.38e-5  3.10e-3  1.00   16262.
10 p[5]        1.08e-3  7.81e-4 0.00101 7.52e-4  7.52e-5  3.10e-3  1.00   15436.
# … with 28 more rows, and 1 more variable: ess_tail <dbl>

Code

mfit$diagnostic_summary()

$num_divergent
[1] 0 0 0 0 0 0 0 0

$num_max_treedepth
[1] 0 0 0 0 0 0 0 0

$ebfmi
[1] 0.9718 0.9143 0.9957 0.9689 0.9794 0.9462 1.0284 0.9763

Trace plots

Code

mcmc_trace(mdrws["beta"])

Code

mcmc_trace(mdrws["alpha"])

Therapeutic dose

Descriptive summary table

save_tex_table(
  make_summary_table(acs_itt_concurc4_dat, "latex"), 
  file.path("outcomes", "secondary", "7-5-anticoagulation-concurrent-therapeutic-summary"))
make_summary_table(acs_itt_concurc4_dat)

Anticoagulation intervention	Patients	Known	Deaths	Any ventilation	DAFV, Median (Q1, Q3)
Low-dose	79	75	3 (4%)	6 (8%)	28 (28, 28)
Intermediate-dose	65	63	1 (2%)	3 (5%)	28 (28, 28)
Therapeutic-dose	50	50	6 (12%)	7 (14%)	28 (28, 28)
Overall	194	188	10 (5%)	16 (8%)	28 (28, 28)

Fit model

ctr <- contr.orthonorm
X <- model.matrix( 
  ~ CAssignment + agegte60 + ctry, 
  data = acs_itt_concurc4_nona_dat, 
  contrasts.arg = list(CAssignment = ctr))[, -1]
y_raw <- ordered(acs_itt_concurc4_nona_dat$out_dafv)
y_mod <- as.integer(y_raw)
N <- dim(X)[1]
K <- dim(X)[2]
unique_y <- length(unique(y_mod))
mdat <- list(
  N = N, K = K, J = unique_y, X = X, y = y_mod, y_raw = y_raw,
  p_par = 2 * rep(1 / unique_y, unique_y),
  beta_sd = c(1, 1, 2.5, 1, 1)
)

snk <- capture.output(
  mfit <- ordmod0$sample(
    data = mdat,
    seed = 49135,
    chains = 8,
    parallel_chains = 8,
    iter_warmup = 1000,
    iter_sampling = 2500,
    refresh = 0
))
mdrws <- as_draws_rvars(mfit$draws())
names(mdrws$beta) <- colnames(X)
mdrws$Ccon <- as.vector(ctr(3) %**% mdrws$beta[1:2])
mdrws$Ctrt <- mdrws$Ccon[-1] - mdrws$Ccon[1]
mdrws$compare <- c("Intermediate vs low" = exp(mdrws$Ctrt[1]),
                   "Therapeutic vs low" = exp(mdrws$Ctrt[2]),
                   "Intermediate vs therapeutic" = exp(mdrws$Ctrt[1] - mdrws$Ctrt[2]))
mdrws$OR <- c("Intermediate" = exp(mdrws$Ctrt[1]),
              "Therapeutic" = exp(mdrws$Ctrt[2]),
              setNames(exp(mdrws$beta[3:5]), c("Age \u2265 60", "India", "Australia/New Zealand")))

Sample distribution

pdat <- acs_itt_concurc4_nona_dat %>%
  dplyr::count(
    CAssignment = factor(
      CAssignment, 
      levels = c("C1", "C2", "C4"),
      labels = c("Low", "Intermediate", "Therapeutic")),
    dafv = ordered(out_dafv, levels = 0:28), 
    .drop = F) %>%
  group_by(CAssignment) %>%
  mutate(p = n / sum(n)) %>%
  mutate(cp = cumsum(p)) %>%
  ungroup()
p <- ggplot(pdat, aes(CAssignment, p)) +
  geom_col(aes(fill = fct_rev(dafv))) +
  labs(x = "Anticoagulation intervention", y = "Cumulative proportion") +
  scale_fill_viridis_d("Days alive and\nfree of ventilation", option = "A", direction = -1) +
  guides(fill = guide_legend(reverse = TRUE, ncol = 3)) +
  theme(legend.key.size = unit(0.75, "lines"))
pth <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(pth, "outcome-7-5-descriptive-concurrent-therapeutic.pdf"), p, height = 2.5, width = 6)
p

Distribution of days alive and free of ventilation by anti-coagulation intervention in subset.

Posterior contrasts

plot_or_densities(mdrws$compare)

Odds ratio summary table

save_tex_table(
  odds_ratio_summary_table_rev(mdrws$OR, "latex"),
  "outcomes/secondary/7-5-primary-model-acs-itt-concurrent-therapeutic-summary-table")
odds_ratio_summary_table_rev(mdrws$OR)

Parameter	Median	95% CrI	Mean (SD)	Pr(OR > 1)
Intermediate	1.46	(0.44, 5.22)	1.81 (1.32)	0.73
Therapeutic	0.57	(0.19, 1.71)	0.67 (0.41)	0.16
Age ≥ 60	0.58	(0.21, 1.61)	0.67 (0.37)	0.15
India	2.63	(0.67, 12.56)	3.64 (3.39)	0.91
Australia/New Zealand	1.13	(0.44, 3.04)	1.28 (0.70)	0.59

Posterior summaries for model parameters.

Diagnostics

Code

mfit$summary()

# A tibble: 24 × 10
   variable       mean   median      sd     mad       q5      q95  rhat ess_bulk
   <chr>         <dbl>    <dbl>   <dbl>   <dbl>    <dbl>    <dbl> <dbl>    <dbl>
 1 lp__       -9.86e+1 -9.83e+1 2.47    2.35    -1.03e+2 -95.2     1.00    8006.
 2 beta_raw[… -6.67e-1 -6.60e-1 0.440   0.432   -1.41e+0   0.0417  1.00   20799.
 3 beta_raw[…  6.76e-2  6.60e-2 0.410   0.408   -6.00e-1   0.747   1.00   21069.
 4 beta_raw[… -2.14e-1 -2.15e-1 0.205   0.205   -5.52e-1   0.122   1.00   13273.
 5 beta_raw[…  9.97e-1  9.66e-1 0.751   0.753   -1.96e-1   2.27    1.00   22147.
 6 beta_raw[…  1.26e-1  1.19e-1 0.494   0.489   -6.73e-1   0.946   1.00   16862.
 7 p[1]        4.45e-2  4.13e-2 0.0193  0.0177   1.90e-2   0.0804  1.00   11411.
 8 p[2]        5.75e-3  4.20e-3 0.00543 0.00390  4.74e-4   0.0165  1.00   18867.
 9 p[3]        5.81e-3  4.22e-3 0.00550 0.00397  4.50e-4   0.0165  1.00   15823.
10 p[4]        5.75e-3  4.14e-3 0.00540 0.00390  4.75e-4   0.0164  1.00   15543.
# … with 14 more rows, and 1 more variable: ess_tail <dbl>

Code

mfit$diagnostic_summary()

$num_divergent
[1] 0 0 0 0 0 0 0 0

$num_max_treedepth
[1] 0 0 0 0 0 0 0 0

$ebfmi
[1] 1.0069 1.0749 0.9556 0.9180 0.9832 0.9993 1.0429 0.9973

Trace plots

Code

mcmc_trace(mdrws["beta"])

Code

mcmc_trace(mdrws["alpha"])